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Abstract 

^^ I We discuss tlie Lebwohl-Laslier model of nematic liquid crystals in a confined geometry, using 

S^ ' Monte Carlo simulation and mean-field theory. A film of material is sandwiched between two 

f— ^ . planar, parallel plates that couple to the adjacent spins via a surface strength e^. We consider the 

LT* I cases where the favoured alignments at the two walls are the same (symmetric cell) or different 



(asymmetric cell). In the latter case, we demonstrate the existence of a single phase transition 
in the slab for all values of the cell thickness. This transition has been observed before in the 
regime of narrow cells, where the two structures involved correspond to different arrangements of 
the nematic director. By studying wider cells, we show that the transition is in fact the usual 
isotropic-to-nematic (capillary) transition under confinement in the case of antagonistic surface 
forces. We show results for a wide range of values of film thickness, and discuss the phenomenology 
using a mean-field model. 

PACS numbers: 61.30.Cz, 61.30.Hn, 61.20. Gy 



I. INTRODUCTION 

The Lebwohl-Lasher lattice spin model [l| is an important model to understand the 
formation of the nematic phase in mesogenic materials. It provides qualitatively correct 
predictions and, in some cases, even quantitative information about nematic properties 
2|, yl . There is renewed interest in the model as regards the behaviour of nematic films and 
the nature of the orientational phase transition. Also recently the confined model has been 
analysed in hybrid geometry [J, |5| . 

In the model, spin unit vectors s are located at the sites of a cubic lattice of lattice 
parameter a. Nearest-neighbour (NN) spins interact via a potential energy — eP2(cos7), 
where e is a coupling parameter (e > 0), and cos 7 = s ■ s', with 7 the relative angle between 
the two spins. P2{x) is the second-degree Legendre polynomial. In the confined model (see 
Fig. [1]), parallel spin layers, h in number, are sandwiched between two planar, parallel plates 
(slit pore geometry), each formed by frozen spins that interact with spins in the first and 
last layers (those adjacent to the plates) also with the same potential, but with a (surface) 
coupling constant e^. The Hamiltonian of the model is then 

■H = -eJ2P2{s-s')-e'^P J2 ^2(^1 ■ ^) - ef ) Yl P^irh, ■ s) (1) 

NN first layer last layer 

where the first sum extends over all distinct NN spins, and the second and third only 
involve the spins in the first and last layers, respectively. The surface coupling constants 
eg , i = 1,2, may or may not be different for both plates. In the simulations to be presented 
below, we take ei and ei to be identical (in Section |V] the case of different constants will 
be considered) but, in general, each plate is assumed to favour a different spin orientation 
(easy axis), rhi or 1712- The case rhi = rh2 is a particular case, the symmetric cell, while 
rhi ^ 1712 is the asymmetric case, also called hybrid or twisted cell, depending on the actual 
orientation of the axes. Since the number of fluctuating spin layers is h, the cell width is 
h + l'vn units of the cubic lattice parameter a. The symmetry of the confined model implies 
that its properties only depend on the scalar mi ■ rh2, and not on the individual components 
of the easy axes. In this respect, our cell is both hybrid (a name reserved for the case where 
one of the axes is normal to its surface, while the other is parallel) and twisted (a situation 
where the two axes lie on the surface planes). 

The situation where rhi ■ it^2 = is very interesting, as the film will be subject to 




FIG. 1: Schematic representation of the confined Lebwohl-Lasher model, h is the number of spin 
layers sandwiched between the two external plates (shaded), /i + 1 is the film thickness in units of 
the cubic lattice parameter. The units vectors along the Cartesian coordinates are indicated. 

antagonistic but equivalent forces at the plates, which create frustration. The nematic 
director can satisfy both surface forces by rotating across the slab, creating an approximately 
linearly dependent, smoothly rotated director configuration (L phase), which involves an 
elastic energy. There have been two recent Monte Carlo (MC) simulations of this model |^,|5|, 
motivated by previous works that indicated the existence of a step-like slab configuration 
(S phase, sometimes called biaxial or exchange-eigenvector phase) in which the director is 
constant except in a thin central region, where it rotates abruptly between the two favoured 
orientations |6|-l9|. These preliminary works, along with a more recent one on the twisted 
cell but with < rhi ■ rh2 < 1 10| , are based on Ginzburg-Landau-type models, and predict 
a L to S (LS) phase transition that was confirmed by the MC studies. A recent analysis of a 
lybrid cell using a surface-force apparatus may have detected this transition experimentally 
ll| . But the nature of the transition, the effect of plate separation, and especially the 
relationship between the LS transition and the bulk behaviour (i.e. isotropic-nematic, or IN 
transition) have not been addressed in MC simulations. Some work on related, continuum 
nematic-fluid slabs under hybrid conditions, analysed by means of density- functional theory, 
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have appeared recently and partially answered some of these questions 

The MC results of Ref . j^ only presented a partial scenario of the problem. As mentioned. 



the connection of the LS transition with the bulk isotropic- nematic phase transition remained 
obscure, and the effect of plate strength e^ and the regime of very small separation were not 
explored. In a recent paper 51], the authors add some confusion to the problem by implicitely 
stating that there is a second transition in the slab, of unknown origin, inferred from a weak 
signal in the specific heat of the slab, as obtained from their MC simulations. 

In the present paper we perform careful MC simulations on the hybrid cell with e^ = e, 
« = 1,2. These simulations will be supplemented by mean-field (MF) theoretical results, 
where cases with ei 7^ ei will also be considered. We obtain the LS phase transition 
from specific-heat data obtained from long MC simulation runs, and extend the analysis to 
very small separations, including the case of a single spin layer. No additional transitions 
are observed in our simulations. The connection with the bulk IN transition is established 
by performing simulations on thicker nematic films, supplemented by MF calculations. The 
available evidence indicates that there is a single transition line in the phase diagram, namely 
the LS transition, and that this transition coincides with the capillary IN transition in the 
confined system, which is connected with the bulk IN transition as the plate separation 

In the remaning sections we first discuss the MC simulation techniques (Section HB . 
and then show the results obtained for the case of symmetric (Section IIIip and asymmetric 
(Section IIVI) plates. The MF model and its results are shown in Section |Vl which includes a 
discussion on the macroscopic approach (Kelvin equation) for this problem. The connection 
with the wetting properties is also discussed. A short discussion on the general picture and 
on relation of the present results with those of Ref. jj] is given in Section IVII Conclusions 
are presented in Section IVIII Some details of the macroscopic model can be found in the 
Appendices. 

II. MONTE CARLO TECHNIQUE 

Let us take each of the h layers to consist oi LxL spins. The total number of spins is then 
A^ = hL"^. The MC simulation runs include two types of moves: one-particle orientational 
moves, and cluster moves. The one-particle orientational moves are carried out using the 

n 

standard algorithm for linear molecules described in Ref. [IJ]. The cluster moves are 



performed by means of the usual bonding criteria for NN particles |l5l . Il6| . The presence of 



the wall-particle interactions imposes some restrictions on the possible reflections that can be 
used to carry out the cluster moves. Notice, however, that the total energy is invariant with 
respect to a simultaneous change of sign of all the x components of the particle orientations. 
The same property applies to the y and z components. Therefore, in our realization of the 
cluster algorithm, we choose at random the component {s^, Sy or Sz) that will eventually 
flip. Then, we test the creation of bonds between every NN pair of particles by taking into 
account the change of interaction energy if only the coordinate of one of the particles of the 
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pair is flipped, the bonding probability being 

hij = 1 - exp <^ min 0, -— (s^jS^jSi • % - s^^^j) \ ', (2) 

where a = {x, y, z} is the chosen direction for the reflections, k is Boltzmann's constant 
and T is the temperature. Once all the possible bonds have been tested, the actual bond 
realization is used to distribute the system in several clusters of particles. The cluster move 
is then performed following the Swendsen-Wang strategy 18| : each cluster is flipped (or not 
flipped) with probability one half. 

The simulations were organized in blocks, each block containing 15000 cycles. A cycle 
consists of trial one-particle orientational moves and one-cluster move. After an equilibration 
period of about 150 blocks, we calculate averages over 175 additional blocks of the potential 
energy per particle u, and the eigenvectors and eigenvalues of different realisations of the 
local Saupe tensor Qi, at each plane i = 1, ...,h. This tensor has components 

iQt)aP = j2 X^ -^i^SakSflk-Sals), a,(] = X,y,Z, (3) 

kGith layer 

where the sum extends over all spins of the ith plane, L^ in number. The local tensor Qi, 
deflned in each layer, is diagonalised, providing eigenvalues Pi, — (Pj— i?j)/2 and — (Pj+i?j)/2. 
The flrst, associated with the x direction in the proper frame (i.e. the frame where Qi is 
diagonal), which coincides with the local nematic director hi, is the local uniaxial nematic 
order parameter, whereas Bi is the biaxial nematic order parameter. The orientation of 
the proper frame with respect to the lab (plate-flxed) frame at each plane is given by the 
tilt angle (pi, which describes the director orientation in the xy plane (spanned by the plate 
orienting fields) and coincides with the angle between the x axes of the two frames. We define 
—IT < (pi < IT. For the symmetric cell, we take rhi = 1712 = x, and (0j) ^ 0. (...) denotes 
a thermal average over spin configurations. For the asymmetric cell, with rhi ■ 1x12 = (we 



take rhi = x and 1x12 = y), we compute, for each layer i, the angle (pi as: 

,2 \ \ 1/2 



COS( 



(<) 



(4) 



where Uxi and Uyi are the x and y components of n^ (thermal averages of local quantities 
at sites lying in the same plane are identical by symmetry). Note that, due to the high 
symmetry of the spin interaction, only one deformation mode of the angle (pi is possible in 
the cell (so that splay, bend and twist are equivalent; see Appendix [U]). To analyze possible 
second-order phase transitions, we also compute an additional order parameter, P^y, with 

N 

y^ SxkSyk ) , (5) 



xy 



k=l 



and, for each plane i, the local order parameters: 



(^= 



xyji 



k^ith layer 



(6) 



The order parameter P^y describes the global orientation of the particles in the plane of the 
interacting fields and is related to the thermal average of the absolute value of one of the 
off-diagonal elements of the Saupe tensor by P^y = (2/3) (IQiyl). Likewise, global uniaxial 
and biaxial order parameters P and B can be defined: 



Notice that these global order parameters do not correspond to those that could be computed 
by diagonalising the global Saupe tensor. The following relation holds between Pi, Bi, {Pxy)i 
and (pi locally (at each plane): 



{P. 



1 



xy/i 



Pi 



El 
3 



|sin20j 



(8) 



Therefore, {Pxy)i reflects the variations of both the nematic order parameters Pi and Bi, 
and of the director tilt angle (pi, across the slab. The LS transition can be monitored in 
principle by the changes with temperature of the global order parameters P, B and P^y As 
we will see, our simulations indicate that the transition in the confined slab has a continuous 
nature in the range of pore widths explored, so that these order parameters do not undergo 
discontinuities, but are singular in their derivatives. The associated singularities are washed 
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out in our (necessarily) finite-size simulations. In fact, the finite-size dependence of the 
order parameters is very weak, and simulations on systems with large lateral sizes, along 
with a proper finite-size scaling analysis, are required. However, relevant response functions 
provide a more clear-cut signature of the transition. We have focused on the excess heat 
capacity per spin, c„ = {du/dT)h, with u = {H) /N the average internal energy per spin. 
In the simulations c^ is obtained from the fluctuations in the energy. The phase-transition 
temperature will be located as that temperature at which c^ reaches a maximum value. 
In order to locate the maximum in the heat capacity we use the synthetic method proposed 



by de Miguel 19|, which we briefly describe in the following. Let us consider that cij are the 
output values of the heat capacity and Acj their associated statistical errors as obtained from 
MC simulations at input temperatures Ti, i = 1, . . .n. Usually we fit cl,- to a polynomial of 
order M in T, Cv{T) = Yli=i o-iT^'^- We search the maximum of this polynomial function by 
computing the value of the temperature, T^, for which the derivative of the heat capacity 
with respect to the temperature is zero, then we calculate c^„ = Cy{Tm). The synthetic 
method consists of the following steps: 

• Generate synthetic sets of n data points, ci- = ci,. + ^j, where ^ is a random number 
drawn from a Gaussian distribution with zero mean value and standard deviation Aq. 

• Find the fitting coefficients a\ and calculate cij corresponding to each synthetic set. 
The set of maximum heat capacities follows a Gaussian distribution, and we determine 
the mean value c^^^{L, h). 

Ik) 

Note that for each synthetic set generated we calculate Tm . This set of temperatures will 
also follow a Gaussian distribution, so we can determine the mean value, which will be 
denoted by Tc,{L, h). 

III. RESULTS FOR THE SYMMETRIC CELL 



First we report on the case of symmetric plates, ifii 



tigated in detail by various authors, using MC simulation 20|, |2l|] , MF theory 2^, l23|] and 



rho. This case has been inves 



renormalisation-group (RG) techniques [2J]. The cases e^ > 0, favouring positive order pa- 
rameter, and e^ < 0, favouring negative order parameter, were considered. Here we focus 
on the first, using ei = ei = e. Mean-field models predict a weak first-order transition, 
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FIG. 2: (Colour online). Excess heat capacity per spin in reduced units, c*, as a function of 
reduced temperature T*, for the symmetric case and for various plate separations (indicated as 
labels). Several lateral sizes, given in the inbox, are considered in each case. 



and a terminal plate separation ht below which the capillary isotropic-nematic transition 
disappears. The plain MF model gives ht = 14, whereas a Bethe model, including two-spin 
correlations |23|, increases the value up to ht = 21. Assuming a monotonic variation due to 
higher-order fluctuations, we may expect ht^ 21. Simulations and RG calculations predict 
a continuous transition, in disagreement with MF results. 

Our own simulation results were based on long runs using the special techniques de- 
scribed in the previous section. Our results, obtained for plate separations h < 24, are not 
compatible with the existence of a phase transition. Fig. [2] shows the behaviour of the heat 
capacity per spin, c* = Cy/k, as a function of reduced temperature T* = kT/e. Various plate 
separations h are shown. In each case an analysis of how the lateral size of the sample L af- 
fects the results has been done. We can see that c^ does not show any significant dependence 
with L (provided that L > /i) as L — t- oo, even for h = 24. Therefore, we may expect ht > 24. 



IV. RESULTS FOR THE HYBRID CELL 

The hybrid cell is the main focus of our work. For this cell we chose rrii = x and m2 = y. 
We have simulated systems with different number of slabs for plate strength ei = ei = e. 
In the following, detailed results are presented for the cases h = 8, which is representative 
of the LS phase transition within a narrow pore, and h = 1, which is a special case. At the 
end of the section the global phase diagram, spanning a wide range of values of h, will be 
discussed. In particular, we show results for the case h = 32, which illustrate the nature 
of the LS phase transition in the regime of wide cells and are used to pinpoint the main 
differences with respect to the regime of narrow cells. 

A. h = 8 

The uniaxial nematic order parameter Pj and the tilt angle (pi profiles are plotted in Fig. 
in] for different values of reduced temperature. In agreement with earlier predictions found 
in the literature [4, |5|, ll2|, ll3|, the orientational structure changes continuously or discontin- 
uously across the slab, depending on the temperature (obviously, one cannot strictly talk 
about continuous or discontinuous functions in a discrete system; these are fuzzy adjectives 
that we ascribe to an interpolating function, passing through all points in the profiles, that 
could reasonably be drawn in each case). For example, the tilt angle clearly shows that, for 
the highest temperature, there is a discontinuity in the centre of the slab, this change be- 
coming steeper as the system size is increased. This is the step-like (S) phase. By contrast, 
at low temperature, the orientation of the director changes smoothly from x to y: this is the 
linear-like (L) phase. At higher or lower temperatures no additional structural changes are 
visible in the order parameters or tilt angle. We conclude that there must be a temperature 
Tc at which the structure changes from the S to the L configuration as a thermodynamic 
phase transition, and that, in view of the smooth variation of the profiles with temperature, 
one can assume this transition to be continuous. Later we will provide evidence that, in the 
thermodynamic limit L — )■ oo, the tilt-angle profile at the transition, corresponding to the 
situation depicted in Fig. M^d), is actually a step function. 

More information about the structural LS transition can be found by looking at the heat 
capacity. The phase transition is signalled by a diverging maximum of the heat capacity as 
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FIG. 3: (Colour online). Neniatic uniaxial order parameter Pi (left panels) and tilt angle of the 
nematic director along the z direction for the slit pore with h = 8. (a) and (b) T* = 0.850; (c) 
and (d) T* = 1.076; (e) and (f) T* = 1.200. The lateral size L used in the simulations is indicated 
in the inbox. 

the system lateral size is increased, Fig. ID^a). The maximum exhibits a linear dependence 
with logL, as shown in Fig. 111(b). This dependence suggests that the confined Lebwohl- 
Lasher system under hybrid conditions for the case h = 8 presents a continuous transition 
belonging to the universality class of the two-dimensional Ising model 25 1. 
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Such a hypothesis is fully supported by considering a cumulant analysis of spin correla- 
tions in the xy plane. Specifically, we define a global Saupe tensor as 

2 = 11:2^' (9) 



and focus on the tensor element Qxy The finite-size dependence 25| of the quantity G4 = 



(2xy)/(2xj/)^ turns out to be fully consistent with the proposed critical behavior. The 
results for h = 8 and different values of L are presented in Fig. O As expected, the different 
curves intersect at values of G4 not too far from the universal value G'4c ~ 1.168 of the 
two-dimensional Ising universality class for systems with L^ = Ly and periodic boundary 



conditions 26|. Then we assume the scaling relation 25 1 

^,(L,/^)=T,(/^) + aL-l/^ (10) 

where the critical exponent has the value u = 1 for the two-dimensional Ising universality 



class (26|, and obtain the critical temperature of the transition as Tc{h) = limL^oo Tc{L, h). 
For the particular value of cell thickness /i = 8 we obtain T*{h = 8) = 1.076(1). 

In Fig. |6]we present results for the quantities Pxy, P, and Pi as a function of temperature 
for the fixed pore width h = 8. In part (a) the dependence of Pxy on lateral size L is shown. 
We see that lateral size hardly affects the value of Pxy in the L phase (low temperatures), 
while the value in the S phase (high temperatures) decreases with lateral size (the location 
of the transition is indicated by an arrow). There is no clear signature of the transition at 
the level of Pxy To check whether Pxy — )■ in the S phase in the thermodynamic limit, 
we have performed extensive simulations for systems with rather large lateral size. The 
results are plotted in Fig. [TJ which represents L^/^Pxy as a function of L~^ (the exponent 
1/8 corresponds to a two-dimensional Ising-like critical transition). From these results one 
can conclude that the transition has a two-dimensional character, at least for the pore size 
h = 8 and smaller (the nature of the transition should change to first order for sufficiently 
wide pores, see discussion in Sec. IIVC|) . 

The uniaxial order parameter. Pi, is plotted in Fig. Et^b) for a fixed lateral size of L = 32 
and for the different planes i = 1,2,3 and 4 (planes with i = 8, 7, 6 and 5 are symmetric). 
The global order parameter P is also plotted. At the transition (indicated by an arrow) the 
order parameter shows a larger variation, but again no anomaly can be seen. Note that the 
variation with temperature is more abrupt for the planes closer to the middle of the pore, 
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FIG. 4: (Colour online), (a) Excess heat capacity per spin in reduced units, c*, for the system with 
/i = 8, as a function of reduced temperature T* and for different lateral system sizes L (indicated 
in the inbox). (b) Maximum of the heat capacity per spin in reduced units, c*'"'"', for the system 
with h = 8, as a function of lateral system sizes L. The straight line is a linear fit. 

which is the region where the director is having more dramatic rearrangements. As opposed 
to Pa;y, the global uniaxial order parameter P should be finite in the thermodynamic limit 
at the transition, and in this limit a kink should exist; again the finite lateral size prevents 
this anomaly to show up. 

The picture that emerges from these results is that, starting from the low-temperature 
region, where the L phase is stable, and on approaching the transition by increasing the 
temperature, the director tilt angle starts to bend from the linear-like configuration and 
ultimately develops an abrupt variation that becomes a step at the transition (so that 
{Pxy)i = at each plane, implying sin20,j = 0). This conclusion is subtle, as it implies that 
the tilt-angle profiles shown in Fig. El^d) for the critical configuration actually tend to a 
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FIG. 5: (Colour online). Dependence of the normalised fourth-order cumulant G4 = {Qty) I {Qlry) 
with reduced temperature T* = kT/e and lateral size L for the system with h = 8. The lateral 
size of the systems is quoted in the legend. 

step-function in the thermodynamic limit L ^> 00. 

The physical nature of the phase transition is easily explained as a competition between 
the anchoring effect of the walls, which the film tries to satisfy simultaneously but creates 
conflicting director orientations at the two walls, and the elastic energy incurred when the 
director rotates between one orientation and the other. At low temperatures or large film 
thickness, the system can accommodate a linearly rotating director in the film. When the 
temperature is high or the film thin, the system prefers to eliminate the (large) elastic 
contribution at the cost of creating a step configuration, which can be regarded as a planar 
defect. 

The L phase is degenerate in the following sense. As one goes from z = Ito z = h through 
a line of sites with equal values of x and y, the orientation of the spins rotates from x to y. 
This rotation can be clockwise (-I-) or anticlockwise (— ). The NN interactions between sites 
impose correlations between pairs of NN site lines, which make favorable that two NN lines 
have the same rotation sign. Below Tc the system chooses (with equal probability) either + 
or — as the preferred orientation sign. 



13 




FIG. 6: (Colour online). Order parameters (a) Pxy and (b) P as a function of temperature T, 
both for pore width h = 8. (a) Different curves give Pxy for different values of lateral size L (see 
key), (b) Order parameter Pi for i = 1 (triangles), i = 2 (open circles), i = 3 (open squares) and 
i = 4 (filled squares) for lateral size L = 32. The global order parameter P is represented by filled 
circles. In both (a) and (b) the vertical arrow indicates the location of the phase transition as 
estimated from the heat capacity. 

B. h=l 

The case h = 1 (single layer) is special. Here the spins are subject to an azimuthally- 
invariant potential that favours spin configurations parallel to the plates. Therefore the 
transition belongs to the XY universality class. In fact, our results for the heat capacity 
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FIG. 7: Dependence of the order parameter P^y on the lateral size L of the system for two different 
values of scaled temperatm'e which are close to the true critical temperature, for the case h = 8. 
Filled circles: T* = 1.076. Open circles: T* = 1.075. Error bars are included in each case. 
The horizontal line indicates an approximate value of L^'^Pxy for the latter temperature in the 
thermodynamic limit L ^>- 00. 

(see Fig. [8], and Table [T]) and the behavior of the nematic order parameter (see Fig. [9]) 
suggest that the transition is of the Berezinskii-Kosterlitz-Thouless (BKT) 27|,|28|] type. The 
heat capacity shows a maximum, but c™^^(L) hardly depends on system size and presents 
a shift towards slightly lower temperatures as L increases. For a given system size L we 
consider the temperature at which \dP/dT\ is maximum as the corresponding pseudo-critical 
temperature Tc{L). With the values for different system sizes a rough estimation of the 
transition temperature in the thermodynamic limit, Tbkt = ^^^^L^ooTdL), can be obtained 
by fitting the results to the equation [29|: 



T,(L)^TBKT + a(logL)- 



With this scheme we obtain T^^^ = T*{h = 1) ~ 0.63 ± 0.01. 



(11) 
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TABLE I: Maximum excess heat capacity per particle, and the corresponding temperatures T^^^, 
and pseudo-critical temperatures, Tc{L), defined as indicated in the text, for the hybrid nematic 
film for pore width h = 1. 



L = 16 L = 24 L = 32 L = 48 L = 64 

c^^''{L)/k 2.661(3) 2.668(3) 2.696(4) 2.628(3) 2.614(3) 

r^^^(L) 0.7210(7) 0.7040(5) 0.6940(7) 0.6872(4) 0.6858(4) 

T*{L) 0.723(1) 0.703(1) 0.691(1) 0.678(1) 0.669(1) 




c;2- 



J"* 



FIG. 8: (Colour online). Variation of the excess heat capacity per particle c* with reduced tem- 
perature T* for h = 1 and different values of L (indicated in the inbox). 

C. Wide pores and global phase diagram 

Using the techniques explained in the previous sections, we have extended the calculation 
of the LS transition to other values of pore width h. The values of Tc{L, h) obtained from 
the heat capacity are used to extrapolate to the thermodynamic limit, using Eqn. ( TTOl) . The 
results of this fitting for the different values of h explored are gathered in Table [III As can 
be seen, the critical temperature Tc{h) increases monotonically with h and approaches the 



value of the bulk isotropic-nematic transition temperature, T{j^ = 1.1225(1) 16|, |30|. One 
important point is that only a single peak is observed in the specific heat in all cases as T 
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FIG. 9: (Colour online). Behavior of the nematic order parameter P as a function of reduced 
temperature T* for h = 1 and various lateral sizes L (indicated in the inbox). 

TABLE II: Estimates of the transition temperatures for the hybrid nematic films. Error bars are 
given between parentheses, in units of the last figure quoted, and correspond to 95% confidence 
level. 



h 
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32 



T*{h) 0.817(1) 0.994(2) 1.076(1) 1.108(1) 1.119(1) 



is varied, indicating the presence of a single transition in this system. Therefore, our data 
do not corroborate the findings of Chiccoh et al. (5|, who claim the existence of two distinct 
peaks in the heat capacity. 

The resulting phase diagram in the plane T-h^^ is presented in Fig. [TOl The interval 
< h~^ < 1 was covered in the MC simulations (the bulk, /i~^ = 0, value was obtained from 



independent simulations in Ref. |16l . |30|). The maximum plate separation considered for 



J] and 



The 



the confined fluid was h = 32, which increases the maximum value used in 
LS transition line spans the whole interval < h~^ < 1. For the plate separations explored, 
1 < h < 32, the transition is continuous. As mentioned before, since the bulk transition is 
of (weakly) first order, there must be a change from first-order to continuous behaviour at 
some (probably large) value of h. 

As the transition line is crossed at fixed h, the spin structure in the slab changes suddenly 
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FIG. 10: Phase diagram for the hybrid cell in the T*-h ^ plane for the case e, 



(1) _ .(2) 



e, showing 



temperatures at which the LS transition occurs for each value of plate separation. L (S) phase is 
stable below (above) the corresponding symbol. Circles: present MC simulation results. Triangle: 
MC simulation results by Chiccoli et al. J]. Squares: present MF results. Filled symbols represent 
first-order phase transitions, while open ones refer to continuous phase transitions. Horizontal 
dotted lines: bulk temperatures as obtained from the MF and MC calculations (upper and lower 
lines, respectively). Continuous line: modified Kelvin equation. Dashed line is a guide to the eye. 



but continuously, as it corresponds to the continuous phase transition discussed in Sec. IIVAI 
Here we show profiles for the cases h = 8 and h = 32, in order to illustrate the differences 
between narrow and wide pores. Fig. [H] shows the change in structure for the cases h = 8 
and 32 as the temperature is increased, reflected by the values of the order parameters Pi 
and {Pxy)i, and by the director tilt angle (pi. At high temperature [Figs. [11] (c) and (f)] the 
structure is of the S type, with an abrupt change in the tilt angle as the middle plane of the 
slab is crossed, and with a low value of the order parameter P in the central region. As T 
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is lowered, we pass from the S to the L structure, with the tilt angle slowly rotating from 
one plate to the other. Note that the value of the order parameter Pi in the slab increases 
substantially at the transition [which occurs in the situations represented by panels (b) and 
(e)]. In the S phase the difference between the two cases shown in the figure, which may be 
representative of a thin {h = 8) and a thick {h = 32) slab, is that, in the thick-slab case, 
the nematic films next to the plates are more separated, leaving a wider orientationally 
disordered region in the central part of the slab. The reason why the central region of the 
pore is not completely disordered, panel (f), may be a finite-size effect. Indeed, as discussed 
in Sec. lIVAI fsee Fig. [7]), we expect a step- function behaviour for (pi at the transition [panel 
(e)] in the thermodynamic limit, while Pi should go to zero right at the middle of the pore, 
and {Pxy)i should be zero everywhere. Therefore, in the situation described in panel (f), the 
tilt-angle profile should be a step function, while the Pi profile would be expected to exhibit 
a wide gap with Pi ^ 0, i.e. a thick isotropic central slab. As h is increased, this central 
region will become thicker, implying that the S phase is the confined phase connected with 
the bulk isotropic phase. On the low-T side, the linear-like phase is the confined nematic 
phase and it evolves to the bulk nematic phase (with a director that rotates more and more 
slowly across the slab). The LS transition is the isotropic-to- nematic (IN) transition in a 
hybrid cell, and no additional capillary or structural transitions should be expected to occur 
in this system. 

As a final comment, we note that the MC data for the transition points obtained by 
Chiccoli et al. j^ are slightly shifted with respect to our own data. These differences 
may result from the more efficient sampling of the present study, which considers cluster 
algorithms in the MC moves. Also, our simulations are longer, maximum lateral sizes are 
larger, and a proper finite-size calculation of the transition temperature is performed. 

V. MEAN-FIELD MODEL 

The MF theor y fo r the Lebwohl-Lasher model has been used before to study symmetric 



nematic slabs 
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L|-|23|. a rich phase diagram with respect to the parameters T, h and surface 
couplings results. Here we use the model to rationalise the MC findings shown in the previous 
section, focusing on the hybrid cell. First we briefiy comment on the implementation of the 
theory and then present the results and their connection with the macroscopic behaviour. 
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FIG. 11: Local order parameters Pi, {Pxy)i and director tilt angle (pi obtained from the MC 
simulations for two different pore widths, h = 8 {L = 64) and h = 32 (L = 48), at various 
temperatures in the neighbourhood of the corresponding critical temperature Tdh). (a) /i = 8 and 
T* = 1; (b) /i = 8 and T* = 1.076; (c) /i = 8 and T* = 1.15; (d) /i = 32 and T* = 1.101; (e) /i = 32 
and T* = 1.118; (f) h = 32 and T* = 1.135. Notice that for T > Tc the jump in (p is system-size 
(L) dependent, and becomes steeper as L approaches the thermodynamic limit. 
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A. Theory and method of solution 

The orientational distribution of a spin in the ith plane is given by the function fi(s). 
The complete MF free-energy functional for the Lebwohl-Lasher model is 

/i „ h 



Fum 



kTj2j dsM^) log [47r/,(s)] - ^ef^ J d^ J ds'Ms)Ms')P2{s ■ s') 



i=l " i=l 



h-1 



- ei'^ J dshis)P,is ■ m,) - ef^ J dsU{s)P2is ■ 7ft,) -J2>^^ J dsf.is). (12) 

The Aj's are Lagrange multipliers ensuring the normalisation J dsfi{s) = 1. The interaction 
part contains contributions from spins on the same layer and from spins on two neighbouring 
layers, and also from the external potentials. Here we find it more convenient to use rhi = x 
and 7712 = i as easy axes. 

Functional minimisation of F provides the corresponding coupled, self-consistent Euler- 
Lagrange equations for each plane, which are projected onto a spherical-harmonics basis 
using 

oo I 

M^) = E E /S^'-(^)- (13) 

/=0 rn=-l 

As usual in MF theory, the corresponding equations can be interpreted as if each spin felt 
an effective field created by their neighbours. The effective field is given by the functions 
(p(")(s), with 

<?(o)(s) = P2{cose), ^W(s) = sin2ecos^, <P^^\s) = sin2ecos2(^, (14) 

and {6, ip) the spherical angles of the spin s. Instead of using the whole distribution functions 
/j(s), the order will be described by three lab-fixed order parameters, rjl where a = 0, 1, 2, 
and i = 1, ...,h runs through the h layers. The order parameters are related to the / = 2- 
subspace coefficients //^ by f^'fj = r]f^ V5/47r, f!^} = -jfli = -r/fV^/evr and f^^ = 
/g _2 = rji yS/Gyr. In terms of 77^ , the Euler-Lagrange equations are written 

r/f) = (^(")(^)>^, ^ = 1,2, ...,/., (15) 
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where (...)^ are averages over the orientational distribution function /j(s), with 



fi{s) oc e "=o 



(16) 



where /3 = 1/kT. In this expression fi{s) has to be normahsed to unity, and we take 

(I 
Vo 

with the properties 



(0) - _i/2, ryW ^ Q^ ^(2) ^ ^^ ^W^ ^ ^^ ^W^ ^ Q ^^^ ^W^ ^ Q_ <?>fH^) are surface fields, 






^(0)(^)_£^(2)(^) 



<?fH^) = <: 



i = h 



1 <i < h. 



(17) 



-(2) 



^(0)(, 



h. 



The order parameters 77^ are related to the eigenvalues of the order tensor, P^ (uniaxial) 
and Bi (biaxial) order parameters, and the director tilt angle (pi, through the relations 



n, 



(0) 



PiP2{cos (j)i) + -Bi sin^ 



^r^=(^r-^ltan: 



1 



r]l'^ = PiSm^(f)i + -Bi{l 



COS^ 0i) 



Here (pi, for the sake of convenience, is measured with respect to the z axis (we remind the 
reader that, due to the symmetry of the model, this angle is the same as the one used in 
the MC simulations). From these equations, we can obtain Pi, Bi and (pi from 7]^ ,77^ and 
?7j- . For the bulk system the surface fields are eliminated and rjl = P, rfi = rfi = 0. 
The isotropic-nematic phase transition is of first order, and occurs at T* = 1.321. The order 
parameter at the transition is P = 0.429. 



B. Results: identical surface couplings 



In this section we consider the confined case and take e^ '' = eY' = e. These values ensure 



,(1) _ .(2) 
-s — ts 

that, at bulk conditions, both surfaces are wet by the nematic phase (see Appendix |X]) so 
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FIG. 12: Order parameters Pi and Sj, and director tilt angle (pi for the two phases coexisting at 
the LS phase transition, in the case e^ = e^ = e and as obtained from MF calculations, (a) L 
phase for h = 8; (b) S phase for h = 8; (c) L phase for h = 9; (d) S phase for h = 9. 

that, close to the bulk transition temperature Ti^, thick nematic films are expected at both 
surfaces. 

Order-parameter and tilt-angle profiles are shown in Fig. [12] for the cases h = 8 and 9 at 
the corresponding transition temperatures. The L and S structures coexist at a first-order 
phase transition, in contrast to the MC results, which indicate a continuous transition. 
As in the case of the MC results deep into the S phase, we note the clear discontinuity 
in the director tilt angle in the coexisting S phase. In the coexisting L phase the director 
configuration adopts a linear-like configuration. Also note that, at the transition, the nematic 
order parameter P changes quite substantially: in the S phase two nematic slabs meet at 
the central region, such that the central spins are almost completely disordered, whereas the 
L phase corresponds to a well-developed nematic slab. The differences between the cases 
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where h is an even or odd number are apparent by comparing the cases h = 8 and h = 9. 
While the L phase hardly changes, the S phase of the even-h case does not have a negligible 
value of the order parameter P at the central region, in contrast with the midpoint of the 
odd-h slab. The biaxial order parameter B is non-negligible only in the neighbourhood of 
the step. The uniaxial order parameter and director tilt angle profiles obtained from the 
MF theory are quite similar to those from MC simulation [cf. the two coexisting phases of 
Figs. [12] (a) and (b) with the structures shown in Figs. [TTT a) and (c)]. 

The nature of the LS transition becomes evident if we look at a wider pore. This is shown 
in Fig. [131 which corresponds to /i = 32. In this case the high-temperature S phase has a 
large central region with a virtually zero value of P. As the transition is crossed from the 
region of high temperature, the value of P in the central region increases to a nematic-like 
value, and the total order parameter in the cell undergoes a discontinuous change. Therefore 
this transition, which corresponds to the capillary isotropic-nematic transition, is the same 
as the LS transition, and it can be concluded that there is a single transition line in the 
phase diagram. Note that the phase corresponding to Fig. [TSTb) (confined isotropic phase, 
with two differently-oriented nematic slabs adsorbed at each plate) for /i = 32 is smoothly 
connected to that of Figs. [T2r b) for h = 8 oi (d) for /i = 9 (the step-like phase), since they 
are actually the same phase but with an 'isotropic' central slab of different width. 

The MF phase diagram, in the plane T vs. h~^, is presented in Fig. [TOl Transition 
temperatures for the different values of h explored are represented by squares. Note that 
the character of the transition changes from first order (for h > 7) to continuous (for h < 6). 
The main differences between the MF results and the MC simulations are: (i) the transition 
is weakly of first order in MF for h > 7; in the simulations it is continuous for the range of 
plate separations explored (as mentioned already, the transition must change to first order 
at some, probably large, value of h, since it is of first order in bulk), (ii) There is a shift 
in the transition to higher values of T in MF, in correspondence with the shift in the bulk 
transition, (iii) In the simulation, the transition line seems to tend to the bulk value from 
below, with a very small slope at the origin h~^ = [see Fig. [10]; in the MF theory, it 
changes slope and actually crosses the bulk temperature at /i ~ 12 (see Fig. [151 where an 
enlarged phase diagram is presented). 

In Fig. [TH we plot the order parameters P and Pxy as a function of reduced temperature 
T* for the case h = 60 [panel (a)], where the LS transition is of first order, and h = 5 [panel 
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FIG. 13: Order parameters Pi and Bi, and director tilt angle (pi for the two phases coexisting for a 
pore width h = 32, as obtained from the MF calculations, in the case es = e^ = e. (a) L phase, 
(b) S phase. 



(b)], where the transition is continuous. The behaviour of the order parameters reflected in 
the flgures may be qualitatively similar to the real situation. In (a), both order parameters 
undergo discontinuous changes, indicated by the dotted vertical lines (the sharp variation 
in P in the metastable step-like branch below the transition corresponds to the frustrated 
wetting transition at each plate due to the confinement). In panel (b) both order parameters 
are continuous but exhibit a 'kink' at the LS transition. Note that Pr^y is always zero in the 
step-like phase, implying that the director tilt-angle is a perfect step function. 
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FIG. 14: Order parameters P and Pxy as a function of reduced temperature T* from mean-field 
theory, (a) h = 60; (b) h = 5. Continuous curves: P. Dashed curves: Pxy In (a), vertical dotted 
lines indicate location of first-order LS phase transition, while arrow points to bulk temperature. 
In (b), arrow indicates temperature of continuous LS transition. 

C. Results: other surface couplings 

We now discuss the case where the surface coupHng constants are different. This situation 
is closer to the experiments. We analyse cases where the couplings of the two surfaces are 
different, and also consider situations where conditions of complete wetting by nematic 
prevail, as well as cases where one or the two surfaces are not wet by the nematic phase. 

• In the first case, the surface couplings are chosen as (ei , ei ) = (e, 0.4e), which again 
ensures a regime of complete wetting by the nematic phase at the two surfaces (see 
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FIG. 15: Phase diagram for the hybrid cell in the t-h~^ plane, with t = {T — Tin)/Tin, for different 
values of the surface couplings e^ and eg (indicated in the inbox), as obtained from mean- field 
theory. Dashed lines correspond to the lowest-order Kelvin equation in each case (see text). 



Appendix |A]). Therefore thick nematic films are expected at both surfaces for temper- 
atures close to the bulk transition temperature. In this case the LS phase-transition 
curve shifts to lower temperatures with respect to the previous case, but by a small 
amount, as evident from Fig. [T3 From a structural point of view, the change involves 
a shift in the location of the step: now it is not symmetrically located with respect 
to the two surfaces, but closer to the surface with the weakest coupling constant (i.e. 
the right surface). This feature can be seen in Figs. fTGT c) and (d), where the uniaxial 
order-parameter profile Pi and director tilt-angle 0j are plotted for the case h = 15 and 
for the two phases coexisting at the LS transition. For comparison, the corresponding 
symmetric profiles for the case ei = ei = e are also plotted in panels (a) and (b). 



(1) A2), 



In the second case, the surface couplings are (e^ , e 



fe, 0.219e). Here conditions 
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FIG. 16: Order parameter Pi and director tilt angle (pi for the two phases coexisting at the LS 
transition for a pore of width h = 15 and different values of the surface coupling constants, as 
obtained from mean-field theory. Upper panels: linear-like phase. Lower panels: step-like phase. 
The surface parameters (eg ,€« ) are as follows: (a) and (b), (e,e); (c) and (d), (e,0.4e); (e) and 
(f), (e,0.219e); and (g) and (h), (0.219e,0.219e). 

of nematic wetting only prevail at one surface (Appendix |A]). The shift in the LS 
transition curve is much more drastic: the maximum is lower, and the curve crosses 
the bulk transition temperature at a higher value of h (see Fig. [T5!) . For narrow pores, 
the profiles now reveal that the step is located next to the weaker surface, as expected 
[see Figs. fTGT e) and (f)]. The pore width at which the transition changes from first to 
second order also moves to higher values (not shown in Fig. [T^ . 

• Finally, we have examined the case (ei , ei ) = (0.219e, 0.219e). Now partial wetting 
applies in both surfaces and, as seen in Fig. [151 the LS transition curve exhibits no 
maximum. The isotropic film in the slab centre is very wide and the two nematic 
films are not in contact except for very narrow pores. The coexistence profiles for the 
linear- and step-like phases for a pore of width h = 15 are plotted in Fig. [T6l(g) and 
(h). 

In summary, as the surface coupling of one of the surfaces is made weaker, the step moves 
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towards that surface, and the LS transition temperature shifts to lower values. The cell width 
where the transition changes from first-order to continuous decreases. The maximum in the 
curve also moves to higher values of pore width, and eventually disappears. This feature 
is related to the wetting properties of the cell, as shown in the macroscopic analysis of the 
following section. 

D. Connection with macroscopic behaviour 

Much of the behaviour shown in the previous sections can be explained using a simple 
macroscopic approach. For a fluid confined into a pore of width /i, the macroscopic Kelvin 
equation gives the undercooling (or overheating) of the transition, with respect to the bulk 
transition, as AT{h) = Tc{h) —Tin = aih~^, with h — )■ oo. As discussed in Appendix [B| the 
coefficient ai can be related to the coexistence parameters A7 and s^ as ai = A^/sj^, where 
Sn is the nematic entropy density at the bulk IN transition, and A7 = 731 — 7sn +7si — 7sn , 
with the superscript denoting the type of substrate, i.e. the left or right substrate (note that 
the value of the surface tensions does not depend on the preferred surface orientation -as 
long as the director remains uniform- but only on the value of the surface coupling e^ ). For 
an isolated surface of type i in contact with a bulk phase, the relation — 7in < 7si — 7sn < 7in 
holds; the right equality corresponds to wetting by nematic, while that in the left pertains 
to wetting by isotropic. Therefore we may have ai < or ai > 0, and the transition curve 
AT{h) will monotonically decrease or increase with h~^, respectively, in the regime of large 
h. The sign of A7 depends on the surface couplings: the difference 7si — 7sn vanishes 
for eg = 0.219e, being positive (negative) for larger (lower) es . For the different surface 
couplings analysed above, we have: 

(i) {ei , es ) = (e, e) and, (e, 0.4e). Since nematic wetting occurs at both surfaces, 731 = 
7sN + 7iN ^see Appendix Rl so that A7 = 271N = 0.0351A;Ta~^. The corresponding ai 
coefficient gives the dashed straight line plotted in Fig. [151 as can be seen, the data 
follow the behaviour predicted by the macroscopic analysis for large h. 

(ii) {el , ei ) = (e, 0.219e). Now wetting occurs only at one surface. Since 731 = 7sn , we 
have A7 = 7iN = 0.0176A;T. Again the MF data follow this behaviour in the regime 
of large h (Fig. [T5|) . 
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(iii) (ei , ei ) = (0.219e, 0.219e). Now partial wetting applies and A7 = 0. The LS tran- 
sition curve departs horizontally from the h axis and therefore exhibits no maximum, 
as indeed shown by the MF results in Fig. [151 

As commented above, the MF results indicate that the two surfaces are wet by the 
nematic phase in the case (ei , ei ) = (e, e). Whether this is also true in the MC simulations 
of Section IIVCI is not known, and a more detailed study of the wetting scenario would be 
necessary. Our present MC data seem to indicate ai < 0, which would be incompatible with 
nematic wetting and even with preferential nematic adsorption, i.e. 7sn < 7si- However, 
since the MC profiles indicate that the plates seem to adsorb preferentially the nematic 
phase, we could still have ai > in the real system, but with a change of regime at very 
large values of h. Another factor to bear in mind is the presumably large correlation length 
of the model, associated with the weakness of the bulk transition, which would give rise 
to slowly decaying interfaces and to the inapplicability of the Kelvin equation except for 
extremely wide pores. 

For smaller separations, the elastic effects in the linear-like phase must be very important, 
since the director rotates essentially between 0° and 90° in a very short distance. These 
effects can be shown (Appendix [B]) to give rise to an additional contribution to the Kelvin 
equation, namely [4] AT(/i) = aih~^ + a2h~'^ , which is the so-called modified Kelvin equation. 
The first term comes from capillary forces, already discussed, whereas the second is due to 
elastic effects. The elastic contribution is always negative (a2 < 0), promoting capillary 
isotropisation, and dominates the physics in the regime of narrow pores. It explains the 
decreasing behaviour of the LS transition curve for narrow pores. The sign of the first term 
dominates for very wide pores, and if positive promotes capillary nematisation, giving rise 
to a maximum in the LS transition curve when combined with the second term. 

In order to estimate the value of 02, it is necessary to compute the elastic constants of 
the model (Appendix [C|) . As shown in Appendix |B| 02 = Kn^/Sst^ = —1.630a^ek~^, with 
K the model elastic constant. Fig. [TOl compares the MF and macroscopic models (note that 
the comparison can only be made in the regime where the transition is of first order). The 
overall agreement is not very good. For large h the surface behaviour correctly predicts the 
capillary LS transition (dashed lines in Fig. [T5|) but, as soon as the pore becomes narrower, 
the elastic contribution comes in. However, because the transition is weakly first-order, the 
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correlation length, and therefore the interfacial thickness, is very large. Consequently, in 
the S phase there are thick nematic layers at the walls, which violate the assumptions of the 
model. In the other cases studied the agreement is also disappointing, in particular in the 
partial-wetting cases. 

VI. DISCUSSION 

The response of the system to confinement is intimately connected to its wetting be- 
haviour. This is especially important in connection with the observation of the step-like 
structure. In a situation of complete wetting of the two surfaces by the nematic phase, the 
thickness of the nematic films at temperatures close to the clearing temperature Tin will 
be large, and the two films with oposing directors will meet at the slab centre when h is 
small, producing a step-like phase which will turn into the linear-like phase as temperature 
is lowered. As the pore gets wider, the step-like phase becomes the isotropic phase with 
a nematic film adsorbed at each surface. In a partial-wetting situation, the nematic film 
thickness will be very small, the central isotropic region will be wide, and the two nematic 
films will never meet, except maybe for very narrow pores. 

In the light of our simulation results and the interpretation obtained from the MF theory, 
it is interesting to discuss the quantity /imax introduced by Chiccoli et al., which these au- 
thors obtain from the intersection between their linearly-extrapolated data for the transition 
temperatures and the bulk temperature Tin (see Fig. 4 in |j]), i.e. the intersection between a 
linear fit to the triangles in our Fig. (TOland the horizontal line T* = 1.1225. Our present MC 
results indicate that the LS transition curve is below the bulk temperature, at least for the 
pore widths explored, and that the transition continues as the confined IN transition up to 
h = oo. Therefore, the value /imax = 16.6 obtained by Chiccoli et al. from the extrapolated 
data, and identified as the maximum slab thickness for which the structural phase transi- 
tion can be found, is somewhat misleading, as it seems to imply that this point terminates 
a phase transition curve; however, the phase transition continues up to /i = oo, the step-like 
phase for narrow pores being smoothly connected (from a thermodynamic viewpoint) with 
the confined isotropic phase for wider pores and eventually with the bulk isotropic phase for 
h = oo. Our results imply that /imax obtained by Chiccoli et al. does not seem to have any 
special meaning |l2l |. 
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VII. CONCLUSIONS 

In this paper we have studied the Lebwohl-Lasher model in a confined sht pore, using 
MC computer simulation and MF theory. Two types of surface conditions have been im- 
posed, namely symmetric and asymmetric walls, with special emphasis on the latter. The 
simulations and the data analysis have been carefully performed, with a view to locating 
accurately the phase transition. For the symmetric walls, we have set a lower limit for the 
pore width at which the capillary isotropic-nematic transition takes place: the transition is 
still absent for h = 24, but the behaviour of the heat capacity indicates that it might occur 
for slightly larger pore widths. 

The asymmetric slab was the central target of our investigations, and consequently was 
studied in more detail. A phase transition, spanning the whole range in pore widths and 
associated with a change from the linear-like to the step-like director configurations (LS 
transition), was measured. In all the cases examined, 1 < h < 32, the transition was 
found to be continuous. The LS transition involves a structural change of the director 
configuration but, from examination of the order-parameter profiles, it is evident that the 
transition corresponds to the IN transition in a confined geometry in a situation where the 
slab is subject to two confiicting favoured directions at the two surfaces. Therefore, there is 
a single phase transition in the confined slab. These results are supported by a MF theory, 
which gives qualitatively similar results. Even though the LS or IN transition is continuous 
in the range 1 < h < 32 according to the simulations, the bulk case, h = oo, presents a 
weakly first-order transition (as obtained from independent simulations). This implies that 
there must be a change in order at some, probably large, value of h; the MF model predicts 
h = 6 but this value is clearly too small. The case h = 1 has also been examined by MC 
simulation. Contrary to the transitions in the case 1 < h < 32, which belong to the 2D 
Ising universality class, when h = 1 the transition is essentially different and pertains to the 
XY-model class. 

The case of different surface coupling constants was also analysed, using only MF theory. 
The results are qualitatively similar, as long as the nematic phase wets both surfaces. In this 
case the IN transition can be more clearly identified with the structural transition studied in 
the literature. When partial wetting applies to one of the surfaces the IN transition occurs 
between two phases, one of which is the linear-like phase; the other, step-like phase, changes 
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in this case to a phase with a director which is uniform in most of the slab volume. When 
neither surface is wet, the latter phase consists of a thick central isotropic slab with thin 
nematic films on the two surfaces. 
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Appendix A: SURFACE TENSIONS AND WETTING PROPERTIES 

The surface tensions of the three interfaces involved are necessary to discuss the wetting 
properties of the model and to investigate how the macroscopic behaviour is obtained in the 
confined slab as /i -> oo. The isotropic-nematic interface 7in was computed in slab geometry, 
by considering a slab of nematic material sandwiched between two isotropic regions at the 
coexistence temperature T* = 1.3212. The uniaxial nematic order-parameter profile is 
depicted in Fig. [T71 The surface tension obtained is 7in = 0.0176kTa~^. From a fit of the 
profile to a hyperbolic-tangent function, we get an interfacial width (correlation length) of 
^ = 2.63a ~ 3a. We note that 7in is relatively small and ^ relatively large, confirming the 
weak character of the bulk isotropic-nematic transition. 

The plate-isotropic, 7si, and plate-nematic, 7sn, surface tensions have also been calculated 
in a range of values of the surface coupling constant e^. Depending on the value of e^, different 
wetting regimes are obtained. For e^ < the wall-nematic has an infinitely thick isotropic 
layer adsorbed on the wall, i.e. 7sn = 7si + 7in, which corresponds to complete wetting of 



the surface- nematic interface by the isotropic phase (cf. simulations results of Ref. (20|). In 
the case e^ > 0.43e the surface tensions satisfy 7si = 7sn + 7in5 implying complete wetting by 
the nematic phase of the surface-isotropic interface. In the interval < e^ < 0.43e a partial 
wetting situation arises. Fig. [18] summarises these results. 
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FIG. 17: Order-parameter profile of the isotropic-neniatic interface in the MF theory for the 
Lebwohl-Lasher model. 

Finally, we note that the sign of the surface-tension difference 7si — 7sn depends on e^. 
It turns out that both 7sn and 7si decrease with e^, but 7si — 7sn < for e^ < 0.219e and 
7si — 7sN > for es > 0.219e. When eg = 0.219e we have 7si — 7sn = 0; this case lies of course 
in the regime of partial wetting (Fig. [T8|) . 

complete wetting partial complete wetting 

by isotropic wetting by nematic 
\ 1 fc. 



t 0.430 
0.219 



8,/8 



FIG. 18: Wetting regime of the Lebwohl-Lasher model in mean-field theory. For e^ < the 
substrate is wet by the isotropic phase. For e^ > 0.430e the substrate is wet by the nematic phase. 
In between a partial- wetting regime occurs. The arrow indicates the case where the surface-nematic 
and surface-isotropic interfaces are equal, which occurs at e^ = 0.219e. 

Appendix B: MACROSCOPIC ANALYSIS 



In order to understand the behaviour of the transition line in the confined system, one 
can use a macroscopic analysis involving capillary and elastic forces and derive a modified 
Kelvin equation. This is valid whenever the transition is of first order. The shift in transition 
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temperature T{h) with respect to the bulk temperature Tin can be obtained by writing 
the free energies of the confined isotropic and nematic phases. For a relative temperature 
AT{h) =T{h) — Ttn but small in absolute value compared to Tj^, we have, for the confined 



isotropic phase [35|: 



5 = 7^'^ + if^ + (//'""' - sAT) h + /., (Bl) 

where //''""'' is the bulk free-energy density of the isotropic phase, Si the entropy density of 
the isotropic phase at IN coexistence, and fs{h) the free energy per unit area of the step 
interface. fs{h) is only appreciable for small /i, i.e. when the two nematic films are in close 
contact, and we neglect it here. 

Now the nematic phase is assumed to consist of a linearly- varying director tilt, with = 
at one wall and = 7r/2 at the other. Then an elastic contribution F^^^^ has to be added: 

^ = 273. + (Z;^--' - .,AT) h+^. (B2) 

At the transition F^ = Fi. Using /.''""'' = //''""^\ and solving for AT: 

.rj. _ 2 (7sN - 7si) , Ktt^ 

where the elastic energy was written as F^i^^ = AhKq'^/2 = AhKir^/Sh'^, with q = 7r/2/i, 
and As = s. — Sj = s. < since Sj = 0. At the transition s. = —OAlSka'^. Using reduced 
units h* = h/a, 7* = ja^/e, T* = kT/e, s* = sa^/k, K* = Ka/e = 3P^ (where P = 0.429 
is the nematic order parameter at the transition, see Table HTB . and 7sn — 7si = — 7in (since 
we are in a wetting situation). 



27;, \ 1 (Sn'P^ 
St ) h* V 8s* J /i*2 



AT* = - Mh + = oMAh*'' - imOh*~\ (B4) 



The first term comes from capillary forces and promotes capillary nematization, whereas the 
second is due to the elastic effects and promotes capillary isotropization. 

Appendix C: ELASTIC CONSTANT 

Since the interaction energy does not couple the relative position of the spins with their 



orientation, there exists no distinction between the three Frank elastic constants 



3l| in 



the Lebwohl-Lasher model, and Ki = K2 = K^ = K. Priest |32| calculated K using a 
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molecular-field theory. Here we rederive the result of Priest in the language of density- 
functional theory, and show numerical results for the bulk elastic constant K. Cleaver and 
Allen 33] have obtained the elastic constant by simulation. 

We consider a smoothy varying director field corresponding to a distorted director. At 
each spin site the director unit vector will point along a different direction, and the excess 
free energy of the bulk, distorted nematic, will be: 

^^^^ = -^*EE I du: jdu,'f{(v.h,)f{u,'.n,)P,{6,.u'). (CI) 

i=l j (NN) 

Here we make explicit the dependence of the distribution function / on the director. Note 
that the director need not be the same on each site. Now consider a smoothy varying director 
field corresponding to a distorted director. We assume that the director rotates about the 
y axis by an angle (f) (see Fig. [1]), with a value proportional to the distance of the jth spin 
from the ith spin along the z axis. Then, assuming <^ 1: 

Cj' ■ hj = w' ■ TZy{(j))ni = {u'^, u'y, u'^) ■ (sin 0, 0, cos 0) = w!. + uj'^(j) - -f4>^ H (C2) 

where hi = (0,0, 1) and 'R-y{(p) is a rotation matrix about the y axis through an angle 0. 
Therefore, 

= fU) + Kf'iu:'.)] 0+1 KV"(u;;) - uj'J'iu',)] 02 + ■ ■ ■ (C3) 

Introducing ( ]C3|) into ( JClj) . subtracting the contribution from the undistorted nematic 
[which is obtained with hi = hj = (0,0,1)], noting that the linear term in vanishes 
by symmetry and that the ideal free-energy term does not contribute to the difference in 
free energy between distorted and undistorted fluid, and going to the continuum by using 
a~^ j dr — !- A^ and — ?- adz(l), one arrives at the following expression for the elastic free 
energy: 

:^^i^ = -^-^jdrjdujdu^'fiuj^) KV"(c.;) - uj'J'iuj'^)] P,{c. ■ (b') [9.0(r)]^(C4) 

z is an effective coordination number, which is the number of neighbours of a given one 
involved in the deformation; since we are rotating about one axis, only 2 out of the 6 
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T* ie*) 


P 


P 


K* 


K* 


K*/P^ 


K*/P^ 




(sim.) 


(theo.) 


(sim.) 


(theo.) 


(sim.) 


(theo.) 


0.400 (2.500) 


0.8922 


0.9260 


2.5290 


2.5726 


3.177 


3.000 


0.750 (1.333) 


0.7672 


0.8406 


1.9747 


2.1200 


3.355 


3.000 


0.900 (1.111) 


0.7668 


0.7901 


1.6448 


1.8730 


3.492 


3.000 


1.000 (1.000) 


0.6863 


0.7471 


1.3103 


1.6745 


3.594 


3.000 


1.080 (0.926) 


0.6038 


0.7041 


0.8587 


1.4872 


3.693 


3.000 


1.321 (0.757) 


— 


0.4290 


— 


0.5520 


— 


3.000 



TABLE III: For different values of scaled temperature T* or inverse scaled temperature e* = e/kT, 
values of uniaxial nematic order parameter P, scaled elastic constant K* , and ratio K* /P"^ from 



33( 1 . The highest temperature corresponds to the bulk phase 



both MF theory and MC simulation 
transition in the MF theory. 

neighbours in the cubic lattice are involved, so that z = 2. To identify the elastic constant, 
we use the expression for the Frank elastic energy 311] • Since our distortion is a bend mode, 
we have, with n = (sin^, 0, cos</)) and cf) = qz (where q is the wavevector of the distortion): 

Fci.. = ^[drK\hx{Vx h)f = \KVq\ (C5) 



where V is the sample volume. Comparing with ( JC4I) . with dz 
expression 



g, we arrive at the 



K* 



K* = Ka/e is the scaled elastic constant 



dd; / duj'f{u) ■ z)P2{u) ■ Lb') {uj' ■ yf f"{uj' ■ z) - {uj' ■ z) f{ib' ■ z) 



(C6) 



?his expression is equivalent to the more general 
one derived by Poniewierski and Stecki 3J] in terms of the direct correlation function. To 
calculate K* , we first expand the distribution function using ( !T3|l and then use the addition 
theorem of spherical harmonics, so that 



\ du}j{(c ■ z)P'2,{ui ■ w') = W— /2oA('^' ■ z). 



(C7) 



Therefore the elastic constant is: 
K* 




[uj' ■ yf fiu' ■ z) - {u' ■ z) f{u' ■ z) P2{u' ■ z). (C8 



37 




-0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 

t 



FIG. 19: Scaled elastic constant K* = Ka/e* for Lebwohl-Lasher model as a function of relative 
temperature t (as defined in caption of Fig. fTOj) . Circles: simulation results of Cleaver and Allen 
20|. Open circles: present MF results. 



Using again the Legendre expansion of the distribution function, taking derivatives, and 
using a couple of recurrence relations for the Legendre polynomials, the integral over a>' can 
be calculated easily, and we obtain the scaled elastic constant: 



K* 



Ait 




P. 



20 



X 




3P' 



(C9) 



where P = {P2{cos6)) is the uniaxial nematic order parameter, P = {P2{cos6)) = 
/2o-\/47r/5. Table IIIII presents a comparison of MF theory with MC simulation. At the 
highest temperature the MF theory overestimates the elastic constant by almost 75 %. The 
temperature T* = 1.08 is 3.8% below the transition temperature from the simulation. At 
the same temperature distance from the MF result (T* = 1.2712) the comparison is quite 
good: K* = 0.8490 from theory versus 0.8587 from simulation. In fact, when plotted versus 
the variable t = {T — Tin)/Tin, the two curves are quite close, see Fig. 
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